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ABSTRACT 

The current understanding of the turbulent dissipation in stellar convective 
zones is based on the assumption that the turbulence follows Kolmogorov scaling. 
This assumption is valid for some cases in which the time frequency of the external 
shear is high (e.g. solar p-modes). However, for many cases of astrophysical 
interest (e.g. binary orbits, stellar pulsations etc.) the timescales of interest 
lie outside the regime of applicability of Kolmogorov scaling. We present direct 
calculations of the dissipation efficiency of the turbulent convective flow in this 
regime, using simulations of anelastic convection with external forcing. We show 
that the effects of the turbulent flow are well represented by an effective viscosity 
coefficient and we provide the values of the effective viscosity as a function of 



the perturbation frequency and co mpare our res u 



for finding the effective viscosity of 



Penevetal 



ts to t he perturbative method 



(l2008bl ) that can be applied to 



actual simulations of the surface convective zones of stars. 



Introduction 



For stars with surface convection, turbulent dissipation in the convective zone is 
believed to be the dominant mechanism responsible for the conversion of mechanical energy 
of tides, stellar oscillations and stellar pulsations to heat. Thus, this is the mechanism 



believ ed to determine the rates of tidal s ynchronization and circularization ( 



1988 



1989) , the amplitudes of st ellar p modes (I Goldreich fc Keeley 



Goldreich et al. 



1977 



Zahn 



1966 



Goldreich fc Kumar 



19941 ). and the instability of stars to pulsations (IGonczi 



19M). 



The simplest approach to estimating the dissipation efficiency of the turbulent 
convection is to assume that the turbulence is homogeneous and isotropic and its interaction 
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with external shear is only local. In that case, we can define an effective viscosity coefficient 
which will capture the effects due to the turbulent flow. The question then is how to find 
the value of this coefficient. 

Currently two prescriptions, based on the assumption of Kolmogorov cascade, exist for 
estimating this coefficient as a function of the time period of the external shear (T). 



Zahnl (11966 



19891 ) proposes that, the effective viscosity should scale linearly with the 



fraction of a turn eddies manage to complete in half a perturbation period: 

1 



v = v max mm 



T 
27 



On the other hand, 



Goldreich fc Nicholson! (119771 ) and 



Goldreich fc Keelevl ( I1977I ). 



argue that eddies with turnover times bigger than T/2ir will not contribute to the 
dissipation. Then Kolmogorov scaling predicts: 



v = v n 



mm 



T 

27TT 



1 



(2) 



Zahn's more efficient dissipation seems to be in better agr eement with observations of 



tidal circularization times for binaries containing a giant star ( Verbunt &: 



the location of the red edge of the Cepheid instability strip (IGonczi 



Phinney 



19951 ) 



19821 ). and even this 



more efficient prescription might be insufficient to ex plain the main sequence circularization 



of binary stars in clusters (IMeibom fc Mathieu 



20051 ). 



However, 



Goldreich fc Keelevl (119771 ) 



Goldreich fc Kumarl ( |l988j ) and 



Goldreich et al 



(119941 ). successfully used the less efficient dissipation to develop a theory for the damping 
of the solar p-modes. In this case the more effective dissipation would require dramatic 
changes in the excitation mechanism in order to explain the observed amplitudes. 



Finally, 



Goodman fc Oh! (119971 ) calculated a lowest order expansion of the effective 



viscosity, which when applied to Kolmogorov turbulence predicts a result closer to the less 
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efficient Goldreich & Nicholson viscosity. While this gives a firmer theoretical foundation 
for the quadratic prescription it does not help with the observational problem of insufficient 
dissipation in the case of tides and stellar pulsations. 

A possible resolution of this problem is suggested by the fact that the successful 
applications of the two prescriptions correspond to very different perturbation periods. The 



Zahnl dl966 



1989) scaling seems to work well for periods of order days, and the Goldreich 



and collaborators quadratic scaling seems to apply for periods of order minutes. This 
distinction is important because, in stars with surface convection, Kolmogorov scaling 
predicts that the eddies with turnover times of several minutes would have typical sizes that 
are very small compared to the local pressure scale height and any other external length 
scales. On the other hand turnover times of days correspond to eddies with typical sizes 
comparable or larger than the local pressure scale height. In this case Kolmogorov scaling 
is not expected to apply. 

The flow seen in 2D and 3D simulations of stellar convection is very different from a 



Kolmogorov cascade (jSofia fe Chan 



1984 



Stein Nordlund 



1989 



Malagoli et al. 



199Q ) 



There are two important distinctions. The first is that the velocity power spectrum is 
much flatter in the simulations than Kolmogorov, and so one expects to find a slower 
loss of dissipation efficiency with increased frequency of the external shear, assuming that 
the dissipation is dominated by the resonant eddies, as long as the external shear has a 
period that corresponds to eddy turnover times too long to fall in the inertial subrange of 
Kolmogorov turbulence. The second is that the flow is no longer isotropic and hence the 
effective viscosity should be a tensor, rather than a scalar. 



As a first at t empt to explore this possibility 



Goodman fe Oh 



Penevetal 



(119971 ) perturbative calculation to the 



(2007. 



Robinson et al. 



2008b) adapted the 



20031 1 numerical 



model of stellar convection and found an asymmetric effective viscosity that scaled linearly 
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with the period of the external shear. However, their perturbative calculation is applicable 
only as long as the forcing period T is small compared to the turnover time of the largest 
eddies. In particular the perturbative treatment is not able to provide the maximum value 
the effective viscosity reaches and the frequency at which it reaches it, which is of great 
importance in calculating tidal interactions and dissipation of pulsations. 



In this article we use the 



PenevetaL 



(j2008al ) spectral anelastic code to perform a direct 



calculation of the turbulent dissipation in a convective zone, by introducing external shear 
as an extra body force in the fluid equations. The goal is to investigate the applicability of 
effective viscosity as an approximation to the actual turbulent dissipa tion, and to deriy e 



direct ly an effective viscosity prescription and compare it against the 



Goodman fc Oh 



( 119971 ) formalism. 



2. Simulations 



2.1. Steady State Convection 



The details of the numerical simulation and the equations evolved are presented in 



PenevetaL 



(j2008al ). We are simulating a rectangular box with impenetrable, constant 
temperature top and bottom boundaries (the z velocity vanishes and the temperature is 
held at some constant value at the top and bottom walls of the box) using the anelastic 
approximation. 

The background state, and the parameters with which all the runs presented in this 
paper were computed , are the same as th e parameters used for all convectively unstable 



tests of section 4.1 of 



PenevetaL 



(l2008al ). For convenience we remind them here: 
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Heat Diffusion Coefficient 

Fig. 1. — The height dependent heat diffusion coefficient 



L x = L y = L z = 4 physical dimensions of the convective box. 

N x = N y = N z = 128 resolution in each direction 

Ptop = 1-0 x 10 5 background pressure at the top of the box 

g = 2.74 gravitational acceleration, in — z direction 

C p = 0.21 specific heat at constant pressure of the fluid 

R = 8.317 x 1CT 2 ideal gas constant of the fluid 

Ti ow = 10.0 temperature at the top boundary of the box 

Thigh — 62.37 temperature at the bottom boundary of the box 

The vertical profile of the heat diffusion coefficient we imposed is presented in figure [TJ 

We initialize the box with random entropy fluctuations and let it evolve with time 
until a steady state is reached. The criteria we used for concluding a steady state has been 
reached were, that the kinetic and thermal energies should stop drifting systematically, and 
only exhibit oscillations at the approximate convective turnover time (see fig. |2J), and that 
the spatial spectra of the velocity and potential temperature remain constant to within a 
few percent. The steady state Fourier spatial and time spectra are presented in fig. [31 
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300 400 

Time 



Fig. 2. — Kinetic (left) and thermal (right) energy content of the convective box used as one 
of the criteria for having reached a steady state. We decided steady state was reached for 
times greater than 400. 

2.2. External Forcing 

After steady state has been reached we introduce external forcing (f) in the form of a 
position and time dependent gravitational acceleration in addition to the already present 
vertical gravity. So, in short, the anelastic momentum equation that we evolve is: 



<9v ~ 9 c 

— = v x u - Vh + -=gz + f , 
at 9 



(3) 



where quantities with tilde represent anelastic perturbations to the background variable 
(denoted by an over bar), v is the velocity vector, h = p/p + v 2 /2 is the enthalpy, 9 is 
the perturbation to the background potential temperature, 9 = f (po/p) R ^ Cp , and f is the 
external forcing. We have investigated the effects of two forms of forcing: 



1. Z (height) dependent horizontal forcing with a Gaussian profile: 



f(z,t) = he 2z cos 



2vrf 



x. 



(4) 




Fig. 3. — The (x: top left, y: top right, z: bottom left, time: bottom right) spectra of the 
3 velocity components and the potential temperature. The thick line in the spatial spectra 
plots corresponds to Kolmogor ov scaling (E^ oc k~ 5 ^ 3 ). The thick line in the time spectra 

(j2008al ) found for the effective viscosity. 



plot corresponds to the scaling 



Penev et al. 
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2. Y dependent forcing in the x direction, again with a Gaussian profile: 



f (y) = he~ 2y2 



cos 



( 



2-Kt 



T 



) 



X. 



(5) 



Clearly the amplitude of the forced velocity will be approximately foT/2iT, so in order 
to investigate the period dependence of the dissipation, we simulated a number of flows with 
different periods and the same foT. That way, the shear created by the external forcing 
was the same for all the flows in the set. We computed two such sets for each forcing case: 
one with f T = 1 and another with f Q T = 0.15. In addition we performed another set of 
simulations with fixed period and varying values of JqT in order to investigate the effects of 
the forcing amplitude on the dissipation efficiency Appendix ?? summarizes the runs and 
the respective number of time steps simulated for each case. 

Notice that the weak forcing cases require a lot more time steps in order to average out 
the turbulent noise and allow us to detect the systematic energy dissipation. Similarly, the 
y dependent forcing requires longer runs as well, because in this case, for the same forcing 
there is less dissipation due to the anisotropic effective viscosity. 

For the strong forcing case, we expect the maximal central velocities to reach 
meixv x = 1/27T ~ 0.16, and for the weak forcing maxv x = 0.15/27T « 0.024. For comparison, 
in our convective box the speed of sound varies between 1.2 at the top of the box and 2.9 at 
the bottom and the typical r.m.s. velocities are between 0.02 and 0.04 (see figure Hj), except 
very near the top and bottom of the box where the collision of the vertical flow with the 
impenetrable boundaries results in higher horizontal velocities. 

In the case of tides in binary stellar systems, the velocities excited by the external 
forcing are small compared to the typical convective velocities, however, performing 
numerical simulations with forcing small enough to ensure that this holds will be prohibitive 
in terms of computational time, because it will require simulations for excessive number 
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0.01 0.02 0.03 0.04 0.05 0.06 0.07 

< V z 2 >!/ 2 



Fig. 4. — The typical steady state r.m.s. velocities (v x top left, v y top right, v z bottom) in 
the absence of forcing. 



of time steps to make sure the dissipation is noticeable among the fluctuations due to the 
turbulence (see section 13721) . 



Results 



3.1. Anisotropic Viscosity 



The viscosity in this work, as in 



PenevetaL 



(j2008af ) is described by a rank-4 tensor 



(Kij mn ), which gives the relation between the strain rate = 1/2(9^- + djVi) and the 
viscous stress a if 



(J; 



■Kijmn&mn i 



(6) 
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where summation over repeated indices is assumed. 

This tensor must obey a set of symmetries in order to avoid infinite torques. In addition 
we expect it to be symmetric with respect of rotations and reflections around the z axis. 
With all these symmetries Kij mn can be shown to have only six independent components. 
Of those, the particular forms of forcing described in section 12721 probe only two families: 

7^1221 = -^2121 = -^2112 =-^1212; 

-^3131 = -^3113 = -^1331 = -^2323 I CO 

/— ^1313) 

-^3232 = -^3223 = -^2332 J 

which will be referred to as X1212 and -^1313 from now on. 



3.2. Direct Viscosity Calculation 

To simplify the discussion we define upfront the following quantities for the z dependent 
forcing: 



C xz {z) = y^v^Xj, yj,z,t w ) cos( 



T 

1,3,111 

,2ttL, 



and S xz (z) = ^ v x (x u yj, z, t w ) sin(— -), (9) 



l,3,iii 



and for the y dependent forcing: 

C xy (y,z) = y^v x (x u y,z,t w )cos(-^-), (10) 

i,w 

and S xy (y,z) = 2J y~i t w )sin(— — ), (11) 

where Xj and yj are the locations of the x and y collocation grid points, and t w are the 
times at which we have sampled the velocity field. We evaluate the sum over an integer 
number of forcing periods T. 
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3.2.1. Depth Dependence of the Effective Viscosity 

We use pre-determined functions of depth, up to a normalization constant, for the 
effective viscosity coefficients that we need, namely: 

#1313 = K° 1313 (v 2 z ) 1/2 H p , (12) 
/ v 2 -\- v 2 \ 

and K 1212 = K° 1212 r^^\ H p . (13) 

This is the same scaling that we used in presenting the perturbative result ?, with the only 
difference that only the velocity in the direction of the external shear is used as the velocity 
scale. In equation [13] we average together both horizontal components of the velocity, 
because on average there should be no physical difference between the two. In practice, the 
different velocity scaling makes little difference, since as we can see from fig. HI away from 
the boundaries all components of the velocity behave alike, except for the fact that v z tends 
to be larger. So, using the full r.m.s. velocity instead of only one component, just leads to 
smaller values of the normalization constants K® 313 and K\ 2V1 . 



3.2.2. Fitting the Spatial Dependence of the Dissipation 

We would like to verify the applicability of the effective viscosity framework to the 
problem of turbulent dissipation, by showing that substituting the turbulent flow with a 
simple viscosity is able to capture not only the total amount of energy dissipated, but also 
the momentum transport, or in other words the spatial dependence of this dissipation. 

For the z dependent x forcing, we would like to show that the work per unit mass done 
by the forcing on the flow at each depth: 

Wf;\z) ee f e- 2z2 C xz (z), (14) 
matches the energy that would be transported and dissipated out of that depth by an 



13 



assumed effective viscosity: 

w:r{z) = Ik 1313 (z)[cuz)c':m+suz)s'L(z)) + 

+ \ (*W*) + ^^1313^)) [C XZ (Z)C' XZ (Z) + S XZ (Z)S' X M) , (15) 

where primes denote derivatives with respect to z and p is the background density profile. 

For the y dependent x forcing, we would like to show that the work per unit mass done 
by the forcing on the flow at each y plane: 

-p Fornax 

W% rb (y) = / p(z)C xy (y, z)dz, (16) 

matches the energy that would be transported and dissipated out of that plane by an 
assumed effective viscosity: 



W™ C (V) = ^ / [C* y (y, z)C%{y, z) + S xy (y, z)S^ y (y, z)] dz, (17) 

where now primes denote derivatives with respect to y, z min and z max are the minimal 
and maximal depth respectively that we want to include in the fit and iV = f Zmax p(z)dz. 
The reason we do not want to include the entire simulated domain is that near the 
boundaries the flow is strongly affected by the impenetrable top and bottom walls and is 
thus non-physical. 

We find the values of A^ 313 and K® 212 by least squares fitting of W x l z sc to W l x u z rh in the 
range z min < z < z max and W X y SC to W X y rb in the range —L y /2 < y < L y /2 respectively. 
Clearly, the presence of the turbulence will cause random fluctuations in the velocity profile 
which should average out if we combine a large enough number of time steps in evaluating 
the quantities C xz , C xy , S xz and S xy . These fluctuations are highly amplified when we 
estimate second derivatives of those quantities, so W x l z sc and W™ y sc suffer much more than 
W x ^ rb and W xy irb . Some representative fitted curves for each series of runs are shown in 
figures El M and [TJ 
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Fig. 5. — Typical least squares fits between W^f b (solid curves) and W™ sc (dashed curves) 
for the two sets of simulations with z dependent x forcing with fixed amplitude and variable 
period. Left: strong forcing, right: weak forcing. In each plot W^ rb and W™ sc have been 
scaled in order to make their maximum values be of order few. 
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Fig. 6. — Typical least squares fits between W^ rb (solid curves) and W™ sc (dashed curves) 
for the two sets of simulations with y dependent x forcing with fixed amplitude and variable 
period. Left: strong forcing, right: weak forcing. In each plot W^ rh and W™ sc have been 
scaled in order to make their maximum values be of order few. 



-16- 



f = 0.8/50 T = 50 f = 0.6/50 T = 50 
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Fig. 7. — The least squares fits between W^ rb (solid curves) and W™ sc (dashed curves) for 
the set of simulations with z dependent x forcing with fixed period and variable amplitude. 
In each plot W^f b and W™ sc have been scaled in order to make their maximum values be of 
order few. 
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We see that, except for one case, the two curves match closely, which shows that the 
effective viscosity assumption captures the effects of turbulent dissipation to a good degree. 
The T = 10 curves for the strong z dependent x forcing case do not match well at all. The 
bad fit is due to the fact that if there is any dissipation present it is undetectable from 
within the turbulent noise. 

The reason for the strongly suppressed dissipation is that for such short periods the 
corresponding eddy sizes are too small to be reliably simulated at our current resolution. 
This fact points to the Goldreich et. al. picture, that the dissipation is dominated by eddies 
with periods close to close that of the external forcing (T), as opposed to the Zahn picture 
in which the largest eddies are always the most important, since in the latter case we would 
expect the linear scaling to continue to arbitrarily short periods. 



3.2.3. Matching Deposited to Dissipated Power 

An alternative way to get a value for the scaling constants K® 313 and K® 212 is to equate 
the overall power deposited into the box by the external forcing: 

/Zmax 
p(z)W t x T\z)dz ,for K° 1313 (18) 
~2"min 

C = N [ LV/2 W%\y)dy , for K° 1212 , (19) 

J-Ly/2 

to the value that the effective viscosity dissipates out of the box: 

p(z)W- c (z)dz for K° 1313 , (20) 
e v J y sc = [ Ly/2 W^(y)dy for K° 1212 . (21) 

J-Ly/2 
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Where we have defined: 



-^1313 (Z 



dC r 



dz 

p(z)K 12 u(z 



dS xz 
dz 

dC X y 

dy 



+ 



dS X y 

dy 



dz. 



(22) 
(23) 



Note that in this case we do not expect to match the spatial dependence, only the 
overall rate. The reason for this, is that viscous forces redistribute the energy in the box as 
well as dissipate it. So we cannot use these quantities to examine the applicability of the 
effective viscosity assumption. 

Because evaluating W^ l z sc (z) and W^y SC (y) requires only first derivatives of the sine 
and cosine velocity components, they suffer significantly less from the turbulent noise 
than W x n z sc {z) and W™ sc (y) from section 13.2.21 On the other hand, some energy transfer 
inevitably occurs near the top and bottom boundaries, which is excluded from evaluating 
the z integrals. For the case of z dependent forcing this is a small amount, since the forcing 
in that part of the box is very small by design. However, for the y dependent forcing 
significant energy transfer does occur in those regions, which could bias the estimated 
viscosities toward lower values (see the last paragraph of section l3~3l) . 

Plots of W^ z sc (z) and W™ c (y) for the same cases as in figures El El and [7] are shown 
in figures [HI M and [TU] respectively, where we have taken z max = —z min = 1.5 for the depth 
dependent forcing and z max = —z min = 1.8 for the y dependent forcing. 



3.3. Comparison Between the Perturbative and Direct Calculation 



We compare the 



Penev et al. 



(l2008al ) perturbative estimate of the ^1313 and K1212 
viscosity coefficients to the directly calculated effective viscosity, obtained by the procedures 
of section I3.2[ in figures [11] and [12j 
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f = 1/300 T = 300 f = 0.15/200 T = 200 




Fig. 8. — The energy dissipation rate due to the effective viscosity at each depth for the 
two sets of simulations with z dependent x forcing with fixed amplitude and variable period. 
Left: strong forcing, right: weak forcing. The plots are for the same simulations as the plots 
in figure El The same scaling has been applied to each plot as in figure [5] 
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f = 1/300 T = 300 f = 0.15/290 T = 290 
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Fig. 9. — The energy dissipation rate due to the effective viscosity at each y plane for the 
two sets of simulations with y dependent x forcing with fixed amplitude and variable period. 
Left: strong forcing, right: weak forcing. The plots are for the same simulations as the plots 
in figure El The same scaling has been applied to each plot as in figure [6] 
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f = 0.8/50 T = 50 f o = 0.6/50 T = 50 
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Fig. 10. — The energy dissipation rate due to the effective viscosity at each depth for the 
set of simulations with z dependent x forcing with fixed period and variable amplitude. The 
plots are for the same simulations as the plots in figure [71 The same scaling has been applied 
to each plot as in figure [7J 
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Fig. 11. — Comparison between the after-the-fact perturbatively estimated i^i 313 (curve) 
and the viscosity obtained from the direct simulations by least squares fit of W™ sc (equation 
M) to W^ rb (equation ED — top ; and from setting f W^ rb (z) = j W^ sc {z) — bottom. 
The strong forcing points are plotted with empty symbols, and the weak forcing with filled 
symbols. The horizontal axis is the perturbation period (T) in units of the convective 
turnover time in the box. Also shown are least squares fits to the strong forcing points 
(dashed line) and the weak forcing points (solid line), by a linear function with saturation. 
The error bars in the top left corner of the plots correspond to the standard deviation of 
K® 313 (assumed the same for all points) obtained from the differences with the fitted curves. 
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Fig. 12. — Comparison between the after-the-fact perturbatively estimated -^1212 (curve) 
and the viscosity obtained from the direct simulations by least squares fit of W™ sc (equation 
ED to W^ rb (equation HE) — top ; and from setting / W%f b (y) = j W™ sc {y) — bottom. 
The strong forcing points are plotted with empty symbols, and the weak forcing with filled 
symbols. The horizontal axis is the perturbation period (T) in units of the convective 
turnover time in the box. Also shown are least squares fits to the strong forcing points 
(dashed line) and the weak forcing points (solid line) by a linear function with saturation. 
The error bars in the top left corner of the plots correspond to the standard deviation of 
K\ 2V1 (assumed the same for all points) obtained from the differences with the fitted curves. 
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The curves in those figures correspond to the perturbative expansion, and the points 
correspond to the direct calculation. Since the perturbative calculation assumes that the 
forcing period is small compared to the convective turnover time, we have used a solid line 
only for T < r c . For longer periods (dotted line) this method is certainly not applicable. 

We see that the effective viscosity predicted with all three methods scales linearly 
with the forcing period for most of the range covered by our simulations. Further, the 
direct calculations show that the effective viscosity saturates for long periods, after which it 
remains roughly constant. One physically expects to see this saturation, because for forcing 
periods much larger than any convective timescale there is no reason why the dissipation 
efficiency of the convective zone should depend on the period. The exact saturation period 
varies between the different forcing cases and even viscosity estimation methods, but it is 
seen to be somewhere in the range 1.5r c < T sat < 3r c . 

In order to obtain a functional dependence for the effective viscosity on period we 
perform a least square fit to the direct calculation points, where we fit a linear function that 
saturates at some period. The resulting fits are shown in figures [TT] and H2] as solid lines for 
the weak forcing cases and dashed lines for the strong forcing cases. The parameters of the 
fitted lines are given in tables [I] and [2j This also allows us to get an estimate of the error 
bars associated with the points. Those are shown in the upper left corners of the plots in 
figures [11] and [12] 

Since for very short periods the perturbative calculation should be valid, and the 
turbulence should be well approximated by a Kolmogorov cascade, we expect that the 
scaling of the effective viscosity with period should be quadratic for those short periods. 
The proportionality constant for this steeper scaling should determine the zero crossings of 
the best fit linear approximations to the effective viscosity. From tables [1] and [2] we see that 
these are mostly positive, except for K® 2 \2 m the weak forcing case, where for the fitting 
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Table 1: The linear regression parameters corresponding to the solid and dashed lines in fig. 



Depth Fit 
strong forcing weak forcing 


Dissipation Match 
strong forcing weak forcing 


slope 


0.10 ±0.012 


0.11 ±0.015 


0.12 ±0.009 


0.11 ±0.02 


zero crossing 


0.086 ±0.013 


0.071 ±0.014 


0.044 ±0.011 


0.02 ±0.02 


saturation period 


2.2 ±0.2 


1.64 ±0.15 


2.4 ±0.2 


2.0 ±0.2 



Table 2: The linear regression parameters corresponding to the solid and dashed lines in fig. 





Depth Fit 


Dissipat 


ion Match 




strong forcing 


weak forcing 


strong forcing 


weak forcing 


slope 


0.097 ±0.023 


0.12 ±0.02 


0.073 ±0.023 


0.11 ±0.01 


zero crossing 


0.065 ±0.023 


-0.003 ±0.023 


0.058 ±0.014 


-0.03 ±0.008 


saturation period 


2.0 ±0.3 


2.2 ±0.2 


2.9 ±0.2 


2.6 ±0.1 



method the zero crossing is basically indistinguishable from zero, and in the dissipation 
matching method it is slightly negative. 

We believe that this negative zero crossing is caused by the fact that the forcing in this 
case is not small near the boundaries, and as a result significant energy is deposited there 
by the external forcing, especially since in those regions the horizontal velocity is very large, 
due to the collision of the mainly vertical flow in the bulk of the box with the impenetrable 
boundaries. This large horizontal flow is non physical and as such we do not include the 
regions near the boundaries in our calculations. Ignoring this external energy source is not 
important for the fitting of the spatial dependence of the deposited energy since it will be 
ignored from both W^ rb and W™ sc . However, if some amount of this energy makes it to 
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the region of the box used to calculate W™ sc before the turbulent cascade has dissipated it, 
it will act to artificially increase W™ sc , and consequently decrease the effective viscosity we 
calculate. 

3.3.1. Amplitude Dependence 

The three dependences in each of the plots of figures dH and [12] correspond to three 
different forcing strengths: strong forcing with Vf orc /v conv ~ 2.7, weak forcing with 
Vforc/vconv ~ 0.4 and the pertubative calculation with Vf orc /v conv <C 1, where Vf orc is the 
peak velocity due to the external forcing and v conv is the root mean square velocity for the 
central plane of the box in the absence of forcing. 

The systematic difference between the three viscosities suggests a possibly important 
amplitude dependence of the dissipation efficiency. In order to get some idea for the 
importance of the magnitude of the external shear in determining the values of K± 313 
we performed four additional simulations with z dependent external forcing with period 
T = r c /2 with strengths intermediate between the strong and weak forcing cases considered 
above (see table ?? for the details of each run and figures [7] and [10] for the energy rate 
curves). Using the same two methods discussed above, we estimate the effective viscosity 
for those cases, which we plot in figure [T3J where we have also added the value of the fitted 
straight lines for the strong and weak forcing from figure [TT] at half the convective turnover 
time. 

Clearly the effective viscosity increases for larger external forcing, but the effect is 
much smaller than the period dependence. Further, the perturbative calculation, which 
would correspond to a forcing strength of zero has a value at T = r c /2 of 0.055, which is far 
below the range of direct calculation values obtained from either method. This is possibly 
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Fig. 13. — The dependence of the effective viscosity K® 313 on the forcing strength 
(max V ext /Vh or ), where Vh or is the horizontal r.m.s. velocity in the absence of forcing. The left 
plot correspond to the effective viscosity estimated by least squares fitting of W™ sc (equation 
[T"5l) to W^ rb (equation [TJj) and the right plot corresponds to the effective viscosity obtained 
by setting / W£ rb (z) = J W™ c (z) 



due to the fact that at the periods covered by our direct calculation points the value of T/r c 
is not small enough for the perturbative calculation to be a good approximation. Another 
possibility is that for weaker forcing the amplitude dependence is much stronger, but we 
must note that according to the perturbative calculation, for very small amplitudes, the 
effective viscosity is independent of the amplitude of the external shear. 

From figure [12] we see that this discrepancy between the direct calculation and the 
perturabative calculation does not exist for Kf 2 u an d the weak forcing direct calculation 
points lie close to the the perturbative calculation curve. 

A possible explanation for the observed amplitude dependence is that, at the forcing 
amplitudes we have used, the external shear is itself going to drive turbulence and lead 
to its own dissipation. Thus, the stronger the forcing, the more additional dissipation is 
expected to occur. 
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4. Conclusions 

We have completed a set of simulations of turbulent convection relevant to stellar 
surface convective zones, with external forcing introduced directly in the momentum 
equation in the form of a periodic, horizontal, position dependent gravitational acceleration. 
We then found an effective viscosity assumed to be of the form of equations [12] and [13] (in 
accordance with mixing length theory), using two methods: 

1. least squares fitting of the work done by the external forcing on the flow at each plane 
of constant forcing to the energy transported or dissipated away from that plane by 
the effective viscosity. 

2. matching the overall energy deposited into the box by the external forcing to the 
energy dissipated by the effective viscosity. 

Both methods produce effective viscosity that scales linearly with period with the same 
slope, but they differ by a constant offset. This offset is due to the fact that, for most 
simulations, the negative side lobes of the energy removed by viscosity are less deep relative 
to the positive central peak than for the external work profile. For the least squares fitting 
method this tends to drive the dimensionless proportionality constant to higher values 
slightly over-predicting the central peak, but getting a bit deeper sidelobes. For the energy 
matching method, this tends to give a lower value for the viscosity in order to match the 
overall integrals under the curves. 

This seems to suggest that the particular depth dependence we assumed for the 
turbulent viscosity might be slightly off. In particular it would seem to require a bit larger 
values at the depths of these side lobes. 

Aside from this small effect, we found that the one parameter fits of the first method 
were able to capture the details of the observed position dependence of the deposited energy 
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(see fig. [8], M and [TO]) , which suggests that the effective viscosity assumption for treating 
turbulent dissipation in convective zones is valid. 



We compared this directly obtained effective 



viscosity wi t 



i the lowest order 



Goodman fc Oh! (119971 ) perturbative expansion of iPenev et al.l (l2008af ). applied to the 
steady state flow without forcing, and we found that this method also predicts linear scaling 
of the effective viscosity with period, having a similar slope to the one predicted by the 
direct calculation. 

Again there is a significant constant offset between the perturbatively calculated 
turbulent viscosity and the viscosity from the direct calculation. Some part of this offset 
might be due to the fact that the perturbative calculation is only valid when the external 
perturbation period is much smaller than the convective turnover time, and the smallest 
value we are able to simulate is T « 0.3r c . 

At least part of this difference is due to actual amplitude dependence of the effective 
viscosity, which can be seen in figure [T3l for K® 313 (t c /2). This amplitude dependence is 
possibly explained by the fact that the external shear caused by the introduced forcing 
is an additional driver of turbulence, and hence dissipation. While our calculations show 
clear evidence of this amplitude dependence, it is relatively less important than the period 
dependence at least in the range of amplitudes accessible to our numerical model. 

The idea that the external shear is acting as an additional driver of turbulence suggest 
an explanation for why we might expect some extra viscosity at the location of the negative 
side lobes of the energy transfer curves. These regions correspond to the locations where 
the shear of the externally forced velocity field is the largest. Thus, if the external forcing 
is driving its own turbulence and causing its own dissipation, we would expect that extra 
viscosity to appear exactly at those locations. However, this explanation would suggest that 
the effect should decrease as the amplitude of the external forcing decreases and we see that 
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this is not the case for the fixed period variable amplitude simulations, where the difference 
between the two methods seems more or less independent of the forcing amplitude. 

To summarize, for forcing periods comparable to the convective turnover time (r c ) the 
effective viscosity scales linearly with period and for our convective box its K1313 and -^1212 
components can be approximated by: 

T 



K 1313 (T) « P (v*) 1/2 H } 



0.23 min 



K 1212 {T) 



P 



2 , 2\ I/ 2 
V x + Vy 



Hn 



2.3r c 
0.23 min 



T 



2.37V. 



1^ + 6 xy (A) 



(24) 
(25) 



where T is the period of the external forcing, and S XZ (A) and 5 xy (A) are constant amplitude 
dependent offsets that are also dependent on the particular method for deriving effective 
viscosities: 





Strong Forcing 




Weak Forcing 




fitted dissipation matching 


fitted 


dissipation matching 




0.10 0.06 


0.08 


0.03 


&xy 


0.09 0.04 


0.02 


-0.02 



The limited spatial resolution of our box does not allow us to reliably simulate the 
case of T r c . In that regime, the assumptions for Kol mogorov turbulence ho ld and the 



effective viscosity should scale quadratically with period ( jGoodman fc Oh 



19971 ). 



We compare the above directly calculated values to a perturbative estimate of the 
effective viscosity, which also predicts linear period dependence with the same slope, but 
with 5 = for all viscosity components. 



Zahnl ( 1966 



19891 ) prescription 



The viscosity of equations [24] and [25] is closer to the 
including the saturation period. However, the fact that we see the dissipation disappearing 
sharply at short periods suggests that the Zahn picture in which the dissipation is dominated 
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by the largest eddies present is not the correct one. If this were the case the linear scaling 
seen at long periods should continue to hold for arbitrarily high frequencies. Instead our 
results suggest that it is the eddies with turnover times does to the external forcing period 
that cause most of the dissipation, and the observed linear scaling is due to the shallower 
than Kolmogorov time spectrum of the convective velocities. As a result, we expect the 
linear loss of efficiency to apply only in a limited range and for shorter periods faster loss of 
efficiency should occur. This might be the resolution of the apparent discrepancy between 
the dissipation necessary to explain tidal circularization and the red edge of the Cepheid 
instability strip on one hand and the observed amplitudes of the solar p-modes on the 
other. For small periods (of order minutes) the Kolmogorov scaling of turbulence holds and 
a quadratic decrease in the effective viscosity is appropriate. For long periods (of order 
days) the assumptions necessary for Kolmogorov cascade are not satisfied and the effective 
viscosity is found to scale linearly with period. 
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